home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / clahe.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  13.9 KB  |  306 lines
  1. /*
  2.  * ANSI C code from the article
  3.  * "Contrast Limited Adaptive Histogram Equalization"
  4.  * by Karel Zuiderveld, karel@cv.ruu.nl
  5.  * in "Graphics Gems IV", Academic Press, 1994
  6.  *
  7.  *
  8.  *  These functions implement Contrast Limited Adaptive Histogram Equalization.
  9.  *  The main routine (CLAHE) expects an input image that is stored contiguously in
  10.  *  memory;  the CLAHE output image overwrites the original input image and has the
  11.  *  same minimum and maximum values (which must be provided by the user).
  12.  *  This implementation assumes that the X- and Y image resolutions are an integer
  13.  *  multiple of the X- and Y sizes of the contextual regions. A check on various other
  14.  *  error conditions is performed.
  15.  *
  16.  *  #define the symbol BYTE_IMAGE to make this implementation suitable for
  17.  *  8-bit images. The maximum number of contextual regions can be redefined
  18.  *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256
  19.  *  contextual regions is not recommended.
  20.  *
  21.  *  The code is ANSI-C and is also C++ compliant.
  22.  *
  23.  *  Author: Karel Zuiderveld, Computer Vision Research Group,
  24.  *         Utrecht, The Netherlands (karel@cv.ruu.nl)
  25.  */
  26.  
  27. #ifdef BYTE_IMAGE
  28. typedef unsigned char kz_pixel_t;     /* for 8 bit-per-pixel images */
  29. #define uiNR_OF_GREY (256)
  30. #else
  31. typedef unsigned short kz_pixel_t;     /* for 12 bit-per-pixel images (default) */
  32. # define uiNR_OF_GREY (4096)
  33. #endif
  34.  
  35. /******** Prototype of CLAHE function. Put this in a separate include file. *****/
  36. int CLAHE(kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes, kz_pixel_t Min,
  37.       kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
  38.       unsigned int uiNrBins, float fCliplimit);
  39.  
  40. /*********************** Local prototypes ************************/
  41. static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
  42. static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,
  43.         unsigned long*, unsigned int, kz_pixel_t*);
  44. static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,
  45.            unsigned int, unsigned long);
  46. static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);
  47. static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,
  48.     unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);
  49.  
  50. /**************     Start of actual code **************/
  51. #include <stdlib.h>             /* To get prototypes of malloc() and free() */
  52.  
  53. const unsigned int uiMAX_REG_X = 16;      /* max. # contextual regions in x-direction */
  54. const unsigned int uiMAX_REG_Y = 16;      /* max. # contextual regions in y-direction */
  55.  
  56.  
  57.  
  58. /************************** main function CLAHE ******************/
  59. int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
  60.      kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
  61.           unsigned int uiNrBins, float fCliplimit)
  62. /*   pImage - Pointer to the input/output image
  63.  *   uiXRes - Image resolution in the X direction
  64.  *   uiYRes - Image resolution in the Y direction
  65.  *   Min - Minimum greyvalue of input image (also becomes minimum of output image)
  66.  *   Max - Maximum greyvalue of input image (also becomes maximum of output image)
  67.  *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
  68.  *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
  69.  *   uiNrBins - Number of greybins for histogram ("dynamic range")
  70.  *   float fCliplimit - Normalized cliplimit (higher values give more contrast)
  71.  * The number of "effective" greylevels in the output image is set by uiNrBins; selecting
  72.  * a small value (eg. 128) speeds up processing and still produce an output image of
  73.  * good quality. The output image will have the same minimum and maximum value as the input
  74.  * image. A clip limit smaller than 1 results in standard (non-contrast limited) AHE.
  75.  */
  76. {
  77.     unsigned int uiX, uiY;          /* counters */
  78.     unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */
  79.     unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */
  80.     unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */
  81.     kz_pixel_t* pImPointer;           /* pointer to image */
  82.     kz_pixel_t aLUT[uiNR_OF_GREY];        /* lookup table used for scaling of input image */
  83.     unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
  84.     unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */
  85.  
  86.     if (uiNrX > uiMAX_REG_X) return -1;       /* # of regions x-direction too large */
  87.     if (uiNrY > uiMAX_REG_Y) return -2;       /* # of regions y-direction too large */
  88.     if (uiXRes % uiNrX) return -3;      /* x-resolution no multiple of uiNrX */
  89.     if (uiYRes & uiNrY) return -4;      /* y-resolution no multiple of uiNrY */
  90.     if (Max >= uiNR_OF_GREY) return -5;       /* maximum too large */
  91.     if (Min >= Max) return -6;          /* minimum equal or larger than maximum */
  92.     if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */
  93.     if (fCliplimit == 1.0) return 0;      /* is OK, immediately returns original image. */
  94.     if (uiNrBins == 0) uiNrBins = 128;      /* default value when not specified */
  95.  
  96.     pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);
  97.     if (pulMapArray == 0) return -8;      /* Not enough memory! (try reducing uiNrBins) */
  98.  
  99.     uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */
  100.     ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;
  101.  
  102.     if(fCliplimit > 0.0) {          /* Calculate actual cliplimit     */
  103.        ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);
  104.        ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;
  105.     }
  106.     else ulClipLimit = 1UL<<14;          /* Large value, do not clip (AHE) */
  107.     MakeLut(aLUT, Min, Max, uiNrBins);      /* Make lookup table for mapping of greyvalues */
  108.     /* Calculate greylevel mappings for each contextual region */
  109.     for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++) {
  110.     for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize) {
  111.         pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];
  112.         MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);
  113.         ClipHistogram(pulHist, uiNrBins, ulClipLimit);
  114.         MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);
  115.     }
  116.     pImPointer += (uiYSize - 1) * uiXRes;          /* skip lines, set pointer */
  117.     }
  118.  
  119.     /* Interpolate greylevel mappings to get CLAHE image */
  120.     for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++) {
  121.     if (uiY == 0) {                      /* special case: top row */
  122.         uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;
  123.     }
  124.     else {
  125.         if (uiY == uiNrY) {                  /* special case: bottom row */
  126.         uiSubY = uiYSize >> 1;    uiYU = uiNrY-1;     uiYB = uiYU;
  127.         }
  128.         else {                      /* default values */
  129.         uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;
  130.         }
  131.     }
  132.     for (uiX = 0; uiX <= uiNrX; uiX++) {
  133.         if (uiX == 0) {                  /* special case: left column */
  134.         uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;
  135.         }
  136.         else {
  137.         if (uiX == uiNrX) {              /* special case: right column */
  138.             uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;
  139.         }
  140.         else {                      /* default values */
  141.             uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 1;
  142.         }
  143.         }
  144.  
  145.         pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];
  146.         pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];
  147.         pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];
  148.         pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];
  149.         Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);
  150.         pImPointer += uiSubX;              /* set pointer on next matrix */
  151.     }
  152.     pImPointer += (uiSubY - 1) * uiXRes;
  153.     }
  154.     free(pulMapArray);                      /* free space for histograms */
  155.     return 0;                          /* return status OK */
  156. }
  157. void ClipHistogram (unsigned long* pulHistogram, unsigned int
  158.              uiNrGreylevels, unsigned long ulClipLimit)
  159. /* This function performs clipping of the histogram and redistribution of bins.
  160.  * The histogram is clipped and the number of excess pixels is counted. Afterwards
  161.  * the excess pixels are equally redistributed across the whole histogram (providing
  162.  * the bin count is smaller than the cliplimit).
  163.  */
  164. {
  165.     unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;
  166.     unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
  167.     long lBinExcess;
  168.  
  169.     ulNrExcess = 0;  pulBinPointer = pulHistogram;
  170.     for (i = 0; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */
  171.     lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;
  172.     if (lBinExcess > 0) ulNrExcess += lBinExcess;      /* excess in current bin */
  173.     };
  174.  
  175.     /* Second part: clip histogram and redistribute excess pixels in each bin */
  176.     ulBinIncr = ulNrExcess / uiNrGreylevels;          /* average binincrement */
  177.     ulUpper =  ulClipLimit - ulBinIncr;     /* Bins larger than ulUpper set to cliplimit */
  178.  
  179.     for (i = 0; i < uiNrGreylevels; i++) {
  180.       if (pulHistogram[i] > ulClipLimit) pulHistogram[i] = ulClipLimit; /* clip bin */
  181.       else {
  182.       if (pulHistogram[i] > ulUpper) {        /* high bin count */
  183.           ulNrExcess -= pulHistogram[i] - ulUpper; pulHistogram[i]=ulClipLimit;
  184.       }
  185.       else {                    /* low bin count */
  186.           ulNrExcess -= ulBinIncr; pulHistogram[i] += ulBinIncr;
  187.       }
  188.        }
  189.     }
  190.  
  191.     while (ulNrExcess) {   /* Redistribute remaining excess  */
  192.     pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;
  193.  
  194.     while (ulNrExcess && pulHisto < pulEndPointer) {
  195.         ulStepSize = uiNrGreylevels / ulNrExcess;
  196.         if (ulStepSize < 1) ulStepSize = 1;          /* stepsize at least 1 */
  197.         for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;
  198.          pulBinPointer += ulStepSize) {
  199.         if (*pulBinPointer < ulClipLimit) {
  200.             (*pulBinPointer)++;     ulNrExcess--;      /* reduce excess */
  201.         }
  202.         }
  203.         pulHisto++;          /* restart redistributing on other bin location */
  204.     }
  205.     }
  206. }
  207. void MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes,
  208.         unsigned int uiSizeX, unsigned int uiSizeY,
  209.         unsigned long* pulHistogram,
  210.         unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)
  211. /* This function classifies the greylevels present in the array image into
  212.  * a greylevel histogram. The pLookupTable specifies the relationship
  213.  * between the greyvalue of the pixel (typically between 0 and 4095) and
  214.  * the corresponding bin in the histogram (usually containing only 128 bins).
  215.  */
  216. {
  217.     kz_pixel_t* pImagePointer;
  218.     unsigned int i;
  219.  
  220.     for (i = 0; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */
  221.  
  222.     for (i = 0; i < uiSizeY; i++) {
  223.     pImagePointer = &pImage[uiSizeX];
  224.     while (pImage < pImagePointer) pulHistogram[pLookupTable[*pImage++]]++;
  225.     pImagePointer += uiXRes;
  226.     pImage = &pImagePointer[-uiSizeX];
  227.     }
  228. }
  229.  
  230. void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,
  231.            unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)
  232. /* This function calculates the equalized lookup table (mapping) by
  233.  * cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].
  234.  */
  235. {
  236.     unsigned int i;  unsigned long ulSum = 0;
  237.     const float fScale = ((float)(Max - Min)) / ulNrOfPixels;
  238.     const unsigned long ulMin = (unsigned long) Min;
  239.  
  240.     for (i = 0; i < uiNrGreylevels; i++) {
  241.     ulSum += pulHistogram[i]; pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);
  242.     if (pulHistogram[i] > Max) pulHistogram[i] = Max;
  243.     }
  244. }
  245.  
  246. void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
  247. /* To speed up histogram clipping, the input image [Min,Max] is scaled down to
  248.  * [0,uiNrBins-1]. This function calculates the LUT.
  249.  */
  250. {
  251.     int i;
  252.     const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);
  253.  
  254.     for (i = Min; i <= Max; i++)  pLUT[i] = (i - Min) / BinSize;
  255. }
  256.  
  257. void Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU,
  258.      unsigned long * pulMapRU, unsigned long * pulMapLB,  unsigned long * pulMapRB,
  259.      unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)
  260. /* pImage      - pointer to input/output image
  261.  * uiXRes      - resolution of image in x-direction
  262.  * pulMap*     - mappings of greylevels from histograms
  263.  * uiXSize     - uiXSize of image submatrix
  264.  * uiYSize     - uiYSize of image submatrix
  265.  * pLUT           - lookup table containing mapping greyvalues to bins
  266.  * This function calculates the new greylevel assignments of pixels within a submatrix
  267.  * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
  268.  * between four different mappings in order to eliminate boundary artifacts.
  269.  * It uses a division; since division is often an expensive operation, I added code to
  270.  * perform a logical shift instead when feasible.
  271.  */
  272. {
  273.     const kz_pixel_t Max = (kz_pixel_t) uiNR_OF_GREY - 1;
  274.     const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
  275.     kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */
  276.  
  277.     unsigned int uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;
  278.  
  279.     if (uiNum & (uiNum - 1))   /* If uiNum is not a power of two, use division */
  280.     for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
  281.      uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
  282.     for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
  283.          uiXCoef++, uiXInvCoef--) {
  284.         GreyValue = pLUT[*pImage];           /* get histogram bin value */
  285.         *pImage++ = (kz_pixel_t ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]
  286.                       + uiXCoef * pulMapRU[GreyValue])
  287.                 + uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
  288.                       + uiXCoef * pulMapRB[GreyValue])) / uiNum);
  289.     }
  290.     }
  291.     else {               /* avoid the division and use a right shift instead */
  292.     while (uiNum >>= 1) uiShift++;           /* Calculate 2log of uiNum */
  293.     for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
  294.          uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
  295.          for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
  296.            uiXCoef++, uiXInvCoef--) {
  297.            GreyValue = pLUT[*pImage];      /* get histogram bin value */
  298.            *pImage++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]
  299.                       + uiXCoef * pulMapRU[GreyValue])
  300.                 + uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
  301.                       + uiXCoef * pulMapRB[GreyValue])) >> uiShift);
  302.         }
  303.     }
  304.     }
  305. }
  306.